logo

Objetivo

Realizar mapas con la mayor granularidad posible, a escala intraurbana, de vulnerabilidad social ante la pandemia de COVID-19. Se pone el foco en la situación de residentes migrantes, grupo para el cual los registros oficiales con frecuencia son escasos o inexistentes.

Con ello se busca contribuir en la toma de decisiones por parte de entidades municipales para la planificación y distribución de atención a nivel intraurbano (víveres con bancos de comida, atención médica in situ, etc.).

En esta instancia se utilizarán solo fuentes abiertas con cobertura global, de modo que el procedimiento sea reproducible en cualquier ciudad de Latinoamérica sin requerir fuentes locales específicas. Por supuesto, los datos locales de gran valor y podrán ser incorporados a éste procedimiento para refinar o extender los resultados.

Esta prueba de concepto toma al asentamiento “Villa Caracas” en Barranquilla, Colombia. Villa Caracas es un asentamiento donde se concentra población migrante. Originalmente un sector de terrenos ociosos, se ha convertido en sitio de edificación de viviendas precarias ante la falta de acceso a vivienda formal.

knitr::include_graphics("img/villa_caracas.png")
Predio del asentamiento Villa Caracas, 2016 (izq.) y 2019 (der.)

Predio del asentamiento Villa Caracas, 2016 (izq.) y 2019 (der.)

Fuentes de datos

Paquetes a utilizar

library(tidyverse) # funciones útiles en general para manipulación y visualización
library(sf) # para datos georeferenciados en formato vectorial
library(osmdata) # para acceso a datos de OpenStreetMap
library(leaflet) # para visualización de datos geográficos
library(nngeo) # para identificar elementos próximos en el espacio  
               # No disponible en CRAN, usar devtools::install_github("michaeldorman/nngeo")

Descarga de datasets

Estimados demográficos

Requerimos estimado georeferenciados de alta resolución espacial.

El portal oficial de datos, para Barranquilla, no ofrece resultados en la categoría “Estadísticas Nacionales”

Tampoco encontramos datos a nivel intra-municipal en el geoportal del Departamento Administrativo Nacional de Estadística (DANE) de Colombia.

Mientras se continúa buscando datos oficiales, podemos recurrir a los estimados demográficos para Colombia producidos por Facebook

url <- "https://data.humdata.org/dataset/2f865527-b7bf-466c-b620-c12b8d07a053/resource/357c91e0-c5fb-4ae2-ad9d-00805e5a075d/download/population_col_2018-10-01.csv.zip"

zipfile <- tempfile()
download.file(url, zipfile)
filename <- unzip(zipfile, list = TRUE)$Name

poblacion <- read_csv(unz(zipfile, filename))

unlink(zipfile)

El dataset contiene puntos definidos por pares de coordenadas en proyección Mercator, y estimados de población en torno a esa posición en 2015 y 2020, para toda Colombia:

head(poblacion)

Cada par latitud/longitud representa un área que corresponde a 1 segundo de arco de resolución (aproximadamente 30m x 30m)

Límites de la ciudad

Pueden descargarse del portal de datos abiertos de Colombia, como “https://www.datos.gov.co/Vivienda-Ciudad-y-Territorio/Limite-Distrito-de-Barranquilla/74fb-mmpu”

limites <- st_read("https://www.datos.gov.co/api/geospatial/74fb-mmpu?method=export&format=GeoJSON")
## Reading layer `OGRGeoJSON' from data source `https://www.datos.gov.co/api/geospatial/74fb-mmpu?method=export&format=GeoJSON' using driver `GeoJSON'
## Simple feature collection with 1 feature and 2 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -74.91979 ymin: 10.91388 xmax: -74.7549 ymax: 11.10681
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
leaflet(limites) %>%
    addProviderTiles(providers$CartoDB.Positron) %>% 
    addPolygons()

Limite del asentamiento Villa Caracas

Debido a su naturaleza reciente e informal, no se espera encontrar un archivo georeferenciado con los límites de Villa Caracas. Como alternativa, se ha trazado y digitalizado su polígono en base a examen de imágenes aérea.

asentamiento <- st_read("data/villa_caracas.geojson")
## Reading layer `villa_caracas' from data source `/home/havb/Documents/repos/vulnerabilidad_COVID/notebook/poc/data/villa_caracas.geojson' using driver `GeoJSON'
## Simple feature collection with 1 feature and 1 field
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -74.81359 ymin: 10.95573 xmax: -74.80969 ymax: 10.96053
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
leaflet(limites) %>%
    addProviderTiles(providers$CartoDB.Positron) %>% 
    addPolygons(fill = NA) %>% 
    addPolygons(data = asentamiento, color = "purple")

Puntos de Interés

Nos interesa estimar accesibilidad de la población a efectores de servicios de salud, como hospitales, centros médicos y farmacias.

En el portal de datos oficial para Barranquilla no se hayan datos en la categoría “Salud y Protección Social”.

Es de esperarse que durante la gestión de la pandemia se identifiquen, o creen, listados oficiales de efectores de salud involucrados en la respuesta sanitarias. Hasta que éstos sean obtenidos, podemos recurrir a información disponible en OpenStreetMap, el repositorio global de información georeferenciada contribuida por voluntarios.

Para obtener datos georeferenciados locales realizaremos consultas a Overpass (http://overpass-api.de/), una interfaz que permite extraer información de la base de datos global de OpenStreetMap.

Overpass requiere que se especifique una “bounding box”, es decir las coordenadas de un rectángulo que abarque la zona de interés (en la práctica, los valores máximos y mínimos de latitud y longitud).

La generamos:

bbox <- matrix(st_bbox(limites), nrow = 2, dimnames = list(c("x", "y"), c("min", "max")))

bbox
##         min       max
## x -74.91979 -74.75490
## y  10.91388  11.10681

Para descargar entidades georeferenciadas se requiere conocer las palabras clave con las que se identifican los registros en la base de OSM. Existe gran detalle para el tipo de datos georeferenciados disponibles: áreas de parques públicos, posición de oficinas de correo o cajeros automáticos, vías de ferrocarril, etc. La nomenclatura se puede consultar en https://wiki.openstreetmap.org/wiki/Map_Features

En este caso vamos a solicitar todas las vías de circulación (calles, avenidas, autopistas, etc) de la ciudad. En la base de datos de OSM ciertas entidades agrupadas bajo la categoría “Healthcare” resultan de inmediato interés:

key value description
amenity clinic A medium-sized medical facility or health centre.
amenity hospital A hospital providing in-patient medical treatment
amenity pharmacy Pharmacy: a shop where a pharmacist sells medications

Para esta prueba de concepto descargaremos clínicas y hospitales.

clinicas <-  opq(bbox) %>% 
  add_osm_feature(key = "amenity", value = "clinic") %>% 
  osmdata_sf() 

hospitales <-  opq(bbox) %>% 
  add_osm_feature(key = "amenity", value = "hospital") %>% 
  osmdata_sf()

Dependiendo del nivel de detalle con el que hayan sido registrados, la mayoría de las clínicas y hospitales están representados por puntos, pero en algunos casos por los polígonos de su planta. Para nuestros fines sólo necesitamos los puntos. Convertiremos los polígonos en puntos, tomando su centroide.

todo_a_puntos <- function(osm_query_sf) {
  

  # extraemos los elementos de la lista que refieren a entidades OSM. Son dataframes con nombre
  # "osm_points", "osm_lines", "osm_polygons", "osm_multilines", "osm_multipolygons"
  puntos <- osm_query_sf[startsWith(names(osm_query_sf), "osm_")] %>% 
    # retiramos los que estan vacíos
    compact() %>% 
    # nos quedamos solo con el campo "name" 
    map(~select(., name)) %>% 
    # combinamos en un solo dataframe
    reduce(rbind) %>% 
    #tomamos los centroides
    st_centroid() 
  
  
  # Solo retenemos los puntos con nombre, ya que los que tienen ese campo vacío
  # parecen ser ubicaciones repetidas
  filter(puntos, !is.na(name))
}
  
hospitales <- cbind(todo_a_puntos(hospitales), clase = "hospital")
clinicas <- cbind(todo_a_puntos(clinicas), clase = "clinica")

hospitales_y_clinicas <- rbind(clinicas, hospitales)

El resultado:

hospitales_y_clinicas
leaflet(hospitales_y_clinicas) %>%
    addProviderTiles(providers$OpenStreetMap) %>% 
    addMarkers(popup = ~paste(clase, name)) %>% 
    addPolygons(data = limites, fill = NA)

Procesamiento

División de la ciudad en zonas de análisis

Aquí es donde cada ciudad puede decidir su unidad geográfica de análisis. Lo ideal serían unidades estadísticas censales, con la mayor resolución (lo más pequeñas) que se disponga.

EL análisis a nivel de agregación similar al barrio (o distritos, comunas, etc) no es deseable debido a su baja resolución espacial. En ausencia de cartografía censal de alta granularidad, puede realizarse la partición de la superficie de la ciudad en celdas arbitrarias, lo suficientemente pequeñas.

Para este ejemplo particionaremos la superficie de la ciudad en celdas de unos 500m de radio (aproximando una superficie de 0.785 km2)

# Usamos una re-proyección equiareal para mayor precisión al calcular superficies

n_celdas <- limites %>% 
  st_transform(crs = "+proj=laea +ellps=WGS84 +units=m +no_defs") %>%
  st_area() %>% 
  {. / (pi * 500^2)} %>% 
  as.numeric() %>% 
  round()

n_celdas
## [1] 199

Definiremos entonces 199 celdas.

celdas <- limites %>%
  st_sample(size = 100000) %>%  # sampleamos al azar un número grande, a gusto
  st_coordinates() %>% # extraemos sendas columnas con long y lat
  kmeans(centers = n_celdas) %>%  # nuestros N grupos de puntos con separación máxima entre si
  .$centers %>% # sólo queremos los centroides
  as_tibble() %>% # convertimos en tibble que es lo que le gusta a sf
  st_as_sf(coords = c("X", "Y"), crs = 4326) %>% # pasamos los centroide a objeto espacial
  st_union() %>% # los combinamos en un sólo objeto multipunto
  st_voronoi() %>% # particionamos en voronoi el espacio donde estan
  st_collection_extract("POLYGON") %>% # del resultado extraemos los polígonos
  st_intersection(limites) %>% # recortamos el cuadrado de los voronoi de acuerdo a los límites
  st_sf %>% # Convertimos en dataframe espacial
  mutate(id = rev(row_number())) # agregamos columna con id
ggplot(celdas) +
  geom_sf() +
  theme_void()

Estimados demogŕaficos por área

Para el ejemplo usaremos los estimados de población mayor a 60 años, como grupo de riesgo. También podría considerase tanto la población general como los mayores, combinando ambas variables un un índice agregado.

En aras de la performance, primero extraemos del dataset nacional los datos que caen dentro de la bounding box de la ciudad (esta operación es muy rápida). Habiendo limitado así los puntos a clasificar, toma mucho menos tiempo verificar a que celda de la ciudad corresponde cada punto del dataset.

poblacion_ciudad <- poblacion %>% 
  filter(between(latitude, bbox["y", "min"], bbox["y", "max"]), 
         between(longitude, bbox["x", "min"], bbox["x", "max"])) %>% 
  st_as_sf(coords = c("longitude", "latitude"), crs = 4326) %>%
  st_join(celdas, left = FALSE) %>% 
  group_by(id)
ggplot(poblacion_ciudad) +
  geom_sf(aes(color = population_2020), size = .0001, alpha = .3) +
  scale_color_viridis_c() +
  theme_minimal() + 
  labs(title = "Población por punto (estimado 2020)", color = "personas")

Agregamos la cantidad de personas por celda, en base a los punto que caen en cada una.

celdas <- poblacion_ciudad %>% 
  st_set_geometry(NULL) %>% 
  group_by(id) %>% 
  summarise(population_2020 = sum(population_2020, na.rm = TRUE)) %>% 
  {left_join(celdas, .)}
ggplot(celdas) +
  geom_sf(aes(fill = population_2020), color = NA) +
  scale_fill_viridis_c() +
  theme_minimal() + 
  labs(title = "Población por área (estimado 2020)", fill = "personas")

Realizamos los mismos pasos para Villa Caracas

poblacion_asentamiento <- poblacion %>% 
  filter(between(latitude, bbox["y", "min"], bbox["y", "max"]), 
         between(longitude, bbox["x", "min"], bbox["x", "max"])) %>% 
  st_as_sf(coords = c("longitude", "latitude"), crs = 4326) %>%
  st_join(asentamiento, left = FALSE)
asentamiento <- poblacion_asentamiento %>% 
  st_set_geometry(NULL) %>% 
  group_by(id) %>% 
  summarise(population_2020 = sum(population_2020, na.rm = TRUE)) %>% 
  {left_join(asentamiento, .)}
asentamiento$population_2020
## [1] 385.0497

Los datos utilizados resultan en un estimado de 385 personas residiendo de Villa Caracas.

Efector de salud más cercano a cada área

Para estimar distancias entre cada celda y su efector de salud más cercano, crearemos una matriz Origen-Destino. Tomamos el centroide de cada celda, y determinamos el efector de salud (hospital o clínica) más cercano a su posición.

Como paso previo, descartamos celdas con población menor a un umbral dado:

umbral <- 1 # Solo celdas con al menos 1 persona

celdas <- celdas %>% 
  filter(population_2020 >= umbral)

Encontramos el elemento de la lista Y (sitios de salud) más cercano a cada elemento de la lista X (centroides de celdas).

En esta instancia tomaremos distancia lineal. En una futura iteración, estimaremos distancia a pie a través de la grilla de calles local usando OSRM

# Calculamos el centroide
centroides <- st_centroid(celdas)

# Identificamos el índice (la posición en el dataframe) de hospital/clínica más cercano a cada uno
id_cercanos <- unlist(st_nn(centroides, hospitales_y_clinicas, progress = FALSE))

# Calculamos la distancia entre cada par (por ahora lineal, a mejorar como distancia a pie por calles)
distancia_a_salud <- st_distance(centroides, hospitales_y_clinicas[id_cercanos,], by_element = TRUE)
# eliminamos la unidad (m)
distancia_a_salud <- as.numeric(distancia_a_salud)

# Agregamos as distancias
celdas <- cbind(celdas, distancia_a_salud)

Y lo mismo para Villa Caracas

# Calculamos el centroide
centroides <- st_centroid(asentamiento)

# Identificamos el índice (la posición en el dataframe) de hospital/clínica más cercano a cada uno
id_cercanos <- unlist(st_nn(centroides, hospitales_y_clinicas))
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |======================================================================| 100%
# Calculamos la distancia entre cada par (por ahora lineal, a mejorar como distancia a pie por calles)
distancia_a_salud <- st_distance(centroides, hospitales_y_clinicas[id_cercanos,], by_element = TRUE)
# eliminamos la unidad (m)
distancia_a_salud <- as.numeric(distancia_a_salud)

# Agregamos as distancias
asentamiento <- cbind(asentamiento, distancia_a_salud) 
ggplot(celdas) +
  geom_sf(aes(fill = distancia_a_salud), color = NA) +
  scale_fill_viridis_c(option = "plasma") +
  theme_minimal() + 
  labs(title = "Distancia a efector de salud más cercano", fill = "m")

Identificación de zonas de riesgo

Combinando los datos demográficos con los de cercanía a salud, resaltamos las zonas donde se concentra población mayor, y a la vez se verifican las mayores distancias.

ggplot(celdas) +
  geom_point(aes(x = population_2020, y = distancia_a_salud), alpha = .5) +
  geom_point(data = asentamiento, aes(x = population_2020, y = distancia_a_salud), color = "purple") +
  labs(x = "personas") +
  theme_minimal()

En el gráfico podemos ver que el punto que corresponde al asentamiento Villa Caracas (en violeta), se encuentra relativamente próximo a un efector de salud.

En base al gráfico de dispersión, consideraremos casos crítico a las celdas con más de 500 residentes mayores que distan más de 1.5 km del centro de salud más cercano

celdas <- celdas %>% 
  mutate(sitio_critico = population_2020 > 500 & distancia_a_salud > 1500)
ggplot(celdas) +
  geom_point(aes(x = population_2020, y = distancia_a_salud)) +
  geom_point(data = filter(celdas, sitio_critico), 
             aes(x = population_2020, y = distancia_a_salud), color = "red", alpha = .5) +
  labs(x = "personas") +
  theme_minimal()

En el mapa:

celdas %>% 
  filter(sitio_critico) %>% 
  leaflet() %>% 
  addProviderTiles(providers$OpenStreetMap) %>% 
  addPolygons(color = "red", popup = ~paste("población:", round(population_2020),"</br>",
                                            "distancia a salud:", round(distancia_a_salud), "m"))